;********************************************WI-CM-1-1-MR-vcbc.nc
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
;load library for  wetbulb_stull
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/heat_stress.ncl"
;The *area_conserve_remap *and *area_conserve_remap_Wrap* have been *deprecated*.
;The* ESMF regrid* <http://www.ncl.ucar.edu/Applications/ESMF.shtml> library is the recommended replacement. It uses a superior methodology.
load "$NCARG_ROOT/lib/ncarg/nclscripts/esmf/ESMF_regridding.ncl"
;********************************************
 begin

;set baseline
 dirPath = "~/source_data/"
 dirPath2 = dirPath+"/outdoor/"

 mo  = (/"ACCESS-CM2","AWI-CM-1-1-MR","BCC-CSM2-MR","CMCC-CM2-SR5","EC-Earth3","GFDL-ESM4", \
             "IITM-ESM","KIOST-ESM","MIROC6","MIROC-ES2L","MPI-ESM1-2-HR","MRI-ESM2-0"/); 12 models
 case = (/"ssp585" /)
 ncase = dimsizes(case);
 nmod = dimsizes(mo)

 syr = 1990; including historical data to identify 1C warming
 fyr = 2099;
 csyr = sprinti("%0.4i",syr);
 cfyr = sprinti("%0.4i",fyr);
 nyr = fyr - syr + 1;
 ndays = nyr*365

;read dimension for 1 deg data
  dimPath=dirPath+"/domain/"
  landfrac = "sftlf_fx_360x180_1deg_gn.nc"
  area = "areacella_fx_360x180_1deg_gn.nc"

  fdim  = addfile(dimPath+landfrac, "r")
  fdim2 = addfile(dimPath+area, "r")
  lat := fdim->lat
  lon := fdim->lon
  nlat = dimsizes(lat);
  nlon = dimsizes(lon);
  garea = fdim2->areacella  ;m^2
  lfrac = fdim->sftlf/100; %
  wgt := garea*1e-6*lfrac ; sqrkm
  wgt@_FillValue = default_fillvalue("float"); make sure FillValue is set correctly for wgt
  copy_VarCoords(garea, wgt)
  printVarSummary(wgt)
  print("global total land area: "+sum(wgt)+" sq km")


  ;--*read population distribution data from "Fig.S2.Population_map_water_rural.ncl"
  ;1) population spending >30mins to access water (rescaled to UNICEF,WHO gross values)
  ;2) rural farmer population data made from 1km SMOD and 1km POP data from EU GHSL
  ;as described in Supplementary Method 2 "Population distribution data"
  dirNc =  dirPath
  filNc2  = "population_data_watercollector_farmer.nc"
  pthNc2  = dirNc+filNc2
  ncdf2 = addfile(pthNc2 ,"r")
  POPR_1deg = ncdf2->POPR_1deg; population engaged in farming
  POPR_dens = ncdf2->POPR_dens
  print("number of agricultural workers in 2020: "+sum(POPR_1deg)*1e-6+" million")

  ;--mask global land area by rural population data
  rural_area = where(POPR_1deg .gt. 0, wgt, wgt@_FillValue); rural area
  print("Rural land area in 2020: "+sum(rural_area)*100*1e-6 +" Mha")


 ;===Use precalculated global warming amount relative to 1980-2009 (based on ERA5, HadCRUT5, BerkeleyEarth, NOAAGlobalTemp) to classify warming levels for each CMIP6 model
  ;see Data_global-warming-amount-ERA5-30yr-CMIP6_12models.ncl
  filename1 = dirPath+"Global_warming_amount_10yrma_"+syr+fyr+"_12models_obs30yr.csv"; 10-year running median
  filename2 = dirPath+"Global_warming_amount_30yrma_"+syr+fyr+"_12models_obs30yr.csv"; 20-year running median

;---Read in file as array of strings so we can parse each line
  lines  = asciiread(filename1,-1,"string")
  lines2 = asciiread(filename2,-1,"string")
  k = 1 ;number of header lines
  nlines = dimsizes(lines)-k   ; skip header lines
  delim = ","
;---Read fields (columns)
  models  =           str_get_field(lines(k:),1,delim)
  exp1   =           str_get_field(lines(k:),2,delim)
  central_10yr =     str_get_field(lines(k:),3,delim)
  period10yr =       str_get_field(lines(k:),5,delim)
 ;read median temperature anomaly from file 1
  Ta_anomaly10yr = tofloat(str_get_field(lines(k:),7,delim))
  ;return to 2d (nmod,nyr); nmod is not 15 but 12 now, be careful
  Ta_anomaly10yr@_FillValue=1e+20

  central_30yr =    str_get_field(lines2(k:),3,delim)
  period30yr =       str_get_field(lines2(k:),5,delim)
  Ta_anomaly30yr = tofloat(str_get_field(lines2(k:),7,delim))
  Ta_anomaly30yr@_FillValue=1e+20



 ;------choose whether 10-yr or 30-yr running period
  Ta_anomaly = Ta_anomaly10yr; either 10-yr or 30-yr running median/max
  periods = period10yr
  central_year = central_10yr

  
  ;=round Ta_anomaly to every 0.1 degC (within an arbitrary tolerance of +-0.05°C)
  ;=to find 10-yr/30-yr median matching given global warming amount most closely (+-0.05°C tolerance, similar to Matthews 2017)
  ;round by 1 decimal to the nearest 0.1 degC (0.951,1.05 to 1, 1.46,1.55 to 1.5, 3.451,3.55 to 3.5, 3.95,4.05 to 4)
  Ta_anomaly01 = decimalPlaces(Ta_anomaly,1,True)
  levs1 = fspan(0.5,6,56) ; 0.1 increment to identify the closest WoE for each grid cell
  nlev1 = dimsizes(levs1)


 ;--warming levels from 0.8 to 3.5 are fully covered by all models using 1990-2099 data
 ;--9/10 models have data for 3.6-4 warming amounts

 ;identify 0.8-4.5 warming amounts for detecting of WoE in the range of 1-4
 levels := fspan(0.8,4.5,38); 0.8 to 4.5 degC by 0.1 degC
 nlev = dimsizes(levels)


 DATA = True; 
 if (DATA) then

 Gday_sun_ens = new((/nmod,nlev,nlat,nlon/),"float")
 Gday_diff_ens = new((/nmod,nlev,nlat,nlon/),"float")
 Gday_diff1_ens = new((/nmod,nlev,nlat,nlon/),"float")
 Gday_diff2_ens = new((/nmod,nlev,nlat,nlon/),"float")
 Gday_zero_ens = new((/nmod,nlev,nlat,nlon/),"float")
 TWday_ens = new((/nmod,nlev,nlat,nlon/),"float")
 Gtwday_ens = new((/nmod,nlev,nlat,nlon/),"float")

 G0_area_sun_ens = new((/nmod,nlev,nlat,nlon/),"float")
 G0_area_diff_ens = new((/nmod,nlev,nlat,nlon/),"float")
 G0_area_zero_ens = new((/nmod,nlev,nlat,nlon/),"float")
 G0_pop_sun_ens = new((/nmod,nlev,nlat,nlon/),"float")
 G0_pop_diff_ens = new((/nmod,nlev,nlat,nlon/),"float")
 G0_pop_zero_ens = new((/nmod,nlev,nlat,nlon/),"float")
 G0_poparea_sun_ens = new((/nmod,nlev,nlat,nlon/),"float")
 G0_poparea_diff_ens = new((/nmod,nlev,nlat,nlon/),"float")
 G0_poparea_zero_ens = new((/nmod,nlev,nlat,nlon/),"float")
 TW35_area_ens = new((/nmod,nlev,nlat,nlon/),"float")
 TW35_pop_ens = new((/nmod,nlev,nlat,nlon/),"float")
 TW35_poparea_ens = new((/nmod,nlev,nlat,nlon/),"float")

 Gday0d_sun_ens = new((/nmod,nlev,nlat,nlon/),"float")
 Gday0d_diff_ens = new((/nmod,nlev,nlat,nlon/),"float")
 Gday0d_zero_ens = new((/nmod,nlev,nlat,nlon/),"float")
 TWday35d_ens = new((/nmod,nlev,nlat,nlon/),"float")
 Gday0d_pop_sun_ens = new((/nmod,nlev,nlat,nlon/),"float")
 Gday0d_pop_diff_ens = new((/nmod,nlev,nlat,nlon/),"float")
 Gday0d_pop_zero_ens = new((/nmod,nlev,nlat,nlon/),"float")
 TWday35d_pop_ens = new((/nmod,nlev,nlat,nlon/),"float")
 G0d_pop_sun_ens = new((/nmod,nlev,nlat,nlon/),"float")
 G0d_pop_diff_ens = new((/nmod,nlev,nlat,nlon/),"float")
 G0d_pop_zero_ens = new((/nmod,nlev,nlat,nlon/),"float")
 TW35d_pop_ens = new((/nmod,nlev,nlat,nlon/),"float")
 RP_lev_ens = new((/nmod,nlev,nlat,nlon/),"float")


 do mm=0,nmod-1

   print(" ")
   print("process CMIP6 model: "+mo(mm))


   infile = dirPath2+"/G-daynight-"+syr+fyr+"-"+mo(mm)+"-vcbc-outdoor.nc"; outdoor work

   fi = addfile(infile, "r")
   Gday_sun  := fi->Gday_mx_sun;  annual maxima of daytime average
   Gday_diff := fi->Gday_mx_diff
   Gday_zero := fi->Gday_mx_zero
   Gday_diff1 := fi->Gday_mx_diff_L
   Gday_diff2 := fi->Gday_mx_diff_U

   Gday0d_sun  := fi->Gday0d_sun
   Gday0d_diff := fi->Gday0d_diff
   Gday0d_zero := fi->Gday0d_zero

   TWday := fi->TWday_mx
   TWday35d := fi->TWday35d


   ;==convert TW to Gtw with constant parameters==
    alpha=0.3  ;reflectance of human body depending on clothes
    rho=1.2 ;density of air (kg/m3)
    cp=1010 ;specific heat capacity of air (Joules/kg/°C)
    pa=101000 ;air pressure (Pa)
    Ts =36
    e_skin = 0.97 ;skin emissivity (Brewster, M. Quinn (1992))
    MH = 300/1.7; moderate metabolic heat outdoors

   ;--use a fixed wind speed to convert TW to Gtw
   ;=Fiala hc function
   u0 =1.9 ; ensemble median U2 at annual maximum of daytime mean TW (not used)
   ; finally dynamic U was used to convert TW to Gtw) in Data.Outdoor_Gday_nomask.ncl
   hc:=14.1*u0^0.5 ;derived from Fig.3 data in Fiala 1999
   ga:=hc/(rho*cp)

   ;Baseline skin temperature for conversion of TW to Gtw
   Ts =35 ;default for TW and Gtw used by Sherwood and Huber (2010)
   esat_skin = satvpr_water_bolton(Ts, (/0,1/))  ;esat at skin temperature 35 degC

   lambda = latent_heat_water(Ts, (/0, 0/), 1, False)
   gamma0 = (pa*cp)/(0.622*lambda)
   esat_tw := satvpr_water_bolton(TWday, (/0,1/))
   ;convert TW to Gtw by Eq.9
   Gtwday := -rho*cp*ga*(Ts - TWday) - rho*cp*ga/gamma0*(esat_skin - esat_tw)
   copy_VarCoords(TWday, Gtwday)


   ;==calc areal and population impacts at each grid cell of each model before summarizing and normalizing by warming level
 
  ;--note it's not medianingful to do ensemble statistics of area and population impacted by Gday>0 at each grid cell, which sets zero area/pop for Gday<0 and exaggerates the ens median or avg area/pop impacts
  ;--whether a grid cell area and population is subject to impact should be determined by ensemble statistics of Gday or TWday

   ;assume the most recent rural popopulation data of 2020 remain constant in future
   RP_3d = conform_dims(dimsizes(Gday0d_sun), POPR_1deg, (/1,2/))
   copy_VarCoords(Gday0d_sun, RP_3d)

   ;2) population-weighted heat stress burden (HSB): person days with Gday>0 at each grid cell
   Gday0d_pop_sun := Gday0d_sun*RP_3d ; person days focus on people who work outdoors on farms 
   
   ;***NOTE to be comparable with outdoor sun, also mask shade/dark scenarios with rural population data
   Gday0d_pop_diff := Gday0d_diff*RP_3d;  
   Gday0d_pop_zero := Gday0d_zero*RP_3d; 
   TWday35d_pop := TWday35d*RP_3d; 


   copy_VarCoords(Gday_sun, Gday0d_pop_sun)
   copy_VarCoords(Gday_sun, Gday0d_pop_diff)
   copy_VarCoords(Gday_sun, Gday0d_pop_zero)
   copy_VarCoords(Gday_sun, TWday35d_pop)
   copy_VarCoords(Gday_sun, RP_3d)


   ;==option 2) 10-yr or 30-yr average 
   ;note Raymond et al. used 99.9th percentile and 30-year maximum, Pal & Elthair and Im et al. used 30-year maximum of daily maximum
   ;we can show 99.9th percentile when comparing with historical period 1980-2019
   ;but use 10-yr/30-yr maximum for each warming level because some models do not have 30-year median Ta for 4degC warming before 2099
 
   ; Set the window size for calc the maximum of a period corresponding to each warming level
   window_size = 10; 20; 30  


   ;==find the 10-yr mean corresponding to each warming level
   do ii=0,nlev-1 ;
     print("")
     print("process warming level: "+levels(ii))
 
     Ta_anomaly_m = Ta_anomaly(ind(models .eq. mo(mm))); select the specific model
     period_m = periods(ind(models .eq. mo(mm))); 10-yr/30-yr periods for each model
     central_year_m = central_year(ind(models .eq. mo(mm)))
 
     tDif = abs(levels(ii) - Ta_anomaly_m)
     ind1 = minind(tDif); closest warming amount
     yr1 = tointeger(central_year_m(ind1))
     ;matching within +-0.05 degC tolerance
     if (.not.(ismissing(ind1)) .and. tDif(ind1) .lt. 0.05) then
       print("closest matching period of "+mo(mm)+" for this level: "+period_m(ind1)+" with "+Ta_anomaly_m(ind1))
       print("at central year: "+yr1); use year index to pick the value when syr of Ta_anomaly is different from model data


      ystart = yr1 - window_size/2 +1
      yend = yr1 + window_size/2
      if (ystart.lt.syr .or. yend.gt.fyr) then
        continue
      end if
      print("start:end year for moving statistics:"+ystart+"-"+yend)

      Gday_sun_ens(mm,ii,:,:) = tofloat(dim_avg_n(Gday_sun({ystart:yend},:,:), 0))
      Gday_diff_ens(mm,ii,:,:) = tofloat(dim_avg_n(Gday_diff({ystart:yend},:,:), 0))
      Gday_diff1_ens(mm,ii,:,:) = tofloat(dim_avg_n(Gday_diff1({ystart:yend},:,:), 0))
      Gday_diff2_ens(mm,ii,:,:) = tofloat(dim_avg_n(Gday_diff2({ystart:yend},:,:), 0))
      Gday_zero_ens(mm,ii,:,:) = tofloat(dim_avg_n(Gday_zero({ystart:yend},:,:), 0))
      Gtwday_ens(mm,ii,:,:) = tofloat(dim_avg_n(Gtwday({ystart:yend},:,:), 0))
      TWday_ens(mm,ii,:,:) = tofloat(dim_avg_n(TWday({ystart:yend},:,:), 0))

      Gday0d_sun_ens(mm,ii,:,:) = tofloat(dim_avg_n(Gday0d_sun({ystart:yend},:,:), 0))
      Gday0d_diff_ens(mm,ii,:,:) = tofloat(dim_avg_n(Gday0d_diff({ystart:yend},:,:), 0))
      Gday0d_zero_ens(mm,ii,:,:) = tofloat(dim_avg_n(Gday0d_zero({ystart:yend},:,:), 0))
      TWday35d_ens(mm,ii,:,:) = tofloat(dim_avg_n(TWday35d({ystart:yend},:,:), 0))
      Gday0d_pop_sun_ens(mm,ii,:,:) = tofloat(dim_avg_n(Gday0d_pop_sun({ystart:yend},:,:), 0))
      Gday0d_pop_diff_ens(mm,ii,:,:) = tofloat(dim_avg_n(Gday0d_pop_diff({ystart:yend},:,:), 0))
      Gday0d_pop_zero_ens(mm,ii,:,:) = tofloat(dim_avg_n(Gday0d_pop_zero({ystart:yend},:,:), 0))
      TWday35d_pop_ens(mm,ii,:,:) = tofloat(dim_avg_n(TWday35d_pop({ystart:yend},:,:), 0))

      ;find the mean population of the 10-yr/30-yr window of each warming level at each grid cell for each model
      ;for later use to calc population impact of the ensemble statistics of Gday>0
      ;Outdoor sun scenario focuses on people who work outdoors 
      RP_lev_ens(mm,ii,:,:) = tofloat(dim_avg_n(RP_3d({ystart:yend},:,:), 0))

     end if
    end do

     ;population with at least one day of Gday>0 (after the 10-year averaging, many values are 0<Gday0d<1 not meaningful)
     G0d_pop_sun_ens(mm,:,:,:) = where(Gday0d_sun_ens(mm,:,:,:).ge.1, RP_lev_ens(mm,:,:,:), RP_lev_ens@_FillValue)
     G0d_pop_diff_ens(mm,:,:,:) = where(Gday0d_diff_ens(mm,:,:,:).ge.1, RP_lev_ens(mm,:,:,:), RP_lev_ens@_FillValue)
     G0d_pop_zero_ens(mm,:,:,:) = where(Gday0d_zero_ens(mm,:,:,:).ge.1, RP_lev_ens(mm,:,:,:), RP_lev_ens@_FillValue)
     TW35d_pop_ens(mm,:,:,:) = where(TWday35d_ens(mm,:,:,:).ge.1, RP_lev_ens(mm,:,:,:), RP_lev_ens@_FillValue)

 end do; nmod



   ;======calculate ensemble median and confidence interval at each grid cell at each warming level
   Gday_sun_med = dim_median_n(Gday_sun_ens, 0)
   Gday_diff_med = dim_median_n(Gday_diff_ens, 0)
   Gday_diff1_med = dim_median_n(Gday_diff1_ens, 0)
   Gday_diff2_med = dim_median_n(Gday_diff2_ens, 0)
   Gday_zero_med = dim_median_n(Gday_zero_ens, 0)
   TWday_med = dim_median_n(TWday_ens, 0)
   Gtwday_med = dim_median_n(Gtwday_ens, 0)

   RP_lev_med = dim_median_n(RP_lev_ens, 0); [nlev,nlat,nlon]

   Gday0d_sun_med = dim_median_n(Gday0d_sun_ens, 0)
   Gday0d_diff_med = dim_median_n(Gday0d_diff_ens, 0)
   Gday0d_zero_med = dim_median_n(Gday0d_zero_ens, 0)
   TWday35d_med = dim_median_n(TWday35d_ens, 0)
   Gday0d_pop_sun_med = dim_median_n(Gday0d_pop_sun_ens, 0)
   Gday0d_pop_diff_med = dim_median_n(Gday0d_pop_diff_ens, 0)
   Gday0d_pop_zero_med = dim_median_n(Gday0d_pop_zero_ens, 0)
   TWday35d_pop_med = dim_median_n(TWday35d_pop_ens, 0)
   G0d_pop_sun_med = dim_median_n(G0d_pop_sun_ens, 0)
   G0d_pop_diff_med = dim_median_n(G0d_pop_diff_ens, 0)
   G0d_pop_zero_med = dim_median_n(G0d_pop_zero_ens, 0)
   TW35d_pop_med = dim_median_n(TW35d_pop_ens, 0)


   ;=========calculate confidence interval at each grid cell 
    print("")
    ncrit = round(nmod/2.0, 3); round to integer 5 or 6 or 7 ;(half of models)
    print("nmod="+nmod+" | ncrit="+ncrit)

   ;--method 2) use the 25-75th percentile interval at each grid cell to account for inter-model uncertainty at specific locations e.g., megacities

   ;show 25–75th percentiles -> 50% confidence interval (often shown in other studies, Zhang et al. 2021, Matthews et al. 2017)
   ;?or show 16-84th percentiles -> 68% confidence interval (similar to 66% the Intergovernmental Panel on Climate Change’s ‘likely range’)

   nval := dim_num_n(.not.ismissing(Gday_sun_ens),0); Counts the number of True models at each warming level
   nval := where(nval.ge.ncrit, nval, default_fillvalue(typeof(nval))); only consider grid cells with >=6 models  
   nval@_FillValue = default_fillvalue(typeof(nval)) 
   printMinMax(nval,1)
   pctl_25 = round(0.25*nval,3)-1; 3rd low -> 3 models (0,1,2) at or below this value
   pctl_75 = round(0.75*nval,3)-1; 9th high -> 9 models (0-8) at or below this value
   ;printVarSummary(pctl_25)
   ;subtract 1 to get the correct model indices

   Gday_sun_sort = Gday_sun_ens
   Gday_diff_sort = Gday_diff_ens
   Gday_diff1_sort = Gday_diff1_ens
   Gday_diff2_sort = Gday_diff2_ens
   Gday_zero_sort = Gday_zero_ens
   TWday_sort = TWday_ens
   Gtwday_sort = Gtwday_ens
   Gday_sun_index = dim_pqsort_n(Gday_sun_sort, 2, 0); sorting models in increasing order at each cell
   Gday_diff_index = dim_pqsort_n(Gday_diff_sort, 2, 0)
   Gday_zero_index = dim_pqsort_n(Gday_zero_sort, 2, 0)
   TWday_index = dim_pqsort_n(TWday_sort, 2, 0)
   Gtwday_index = dim_pqsort_n(Gtwday_sort, 2, 0)

   Gday0d_sun_sort = Gday0d_sun_ens
   Gday0d_diff_sort = Gday0d_diff_ens
   Gday0d_zero_sort = Gday0d_zero_ens
   TWday35d_sort = TWday35d_ens
   Gday0d_sun_index = dim_pqsort_n(Gday0d_sun_sort, 2, 0); sorting models in increasing order at each cell
   Gday0d_diff_index = dim_pqsort_n(Gday0d_diff_sort, 2, 0)
   Gday0d_zero_index = dim_pqsort_n(Gday0d_zero_sort, 2, 0)
   TWday35d_index = dim_pqsort_n(TWday35d_sort, 2, 0)
   Gday0d_pop_sun_sort = Gday0d_pop_sun_ens
   Gday0d_pop_diff_sort = Gday0d_pop_diff_ens
   Gday0d_pop_zero_sort = Gday0d_pop_zero_ens
   TWday35d_pop_sort = TWday35d_pop_ens
   Gday0d_pop_sun_index = dim_pqsort_n(Gday0d_pop_sun_sort, 2, 0); sorting models in increasing order at each cell
   Gday0d_pop_diff_index = dim_pqsort_n(Gday0d_pop_diff_sort, 2, 0)
   Gday0d_pop_zero_index = dim_pqsort_n(Gday0d_pop_zero_sort, 2, 0)
   TWday35d_pop_index = dim_pqsort_n(TWday35d_pop_sort, 2, 0)


   Gday_sun_low = new((/nlev,nlat,nlon/), typeof(Gday_sun))
   Gday_diff_low = new((/nlev,nlat,nlon/), typeof(Gday_sun))
   Gday_zero_low = new((/nlev,nlat,nlon/), typeof(Gday_sun))
   TWday_low = new((/nlev,nlat,nlon/), typeof(TWday))
   Gtwday_low = new((/nlev,nlat,nlon/), typeof(Gtwday))
   Gday_diff_high = new((/nlev,nlat,nlon/), typeof(Gday_sun))
   Gday_sun_high = new((/nlev,nlat,nlon/), typeof(Gday_sun))
   Gday_zero_high = new((/nlev,nlat,nlon/), typeof(Gday_sun))
   TWday_high = new((/nlev,nlat,nlon/), typeof(TWday))
   Gtwday_high = new((/nlev,nlat,nlon/), typeof(Gtwday))

   Gday0d_sun_low = new((/nlev,nlat,nlon/), typeof(Gday0d_sun))
   Gday0d_diff_low = new((/nlev,nlat,nlon/), typeof(Gday0d_sun))
   Gday0d_zero_low = new((/nlev,nlat,nlon/), typeof(Gday0d_sun))
   TWday35d_low = new((/nlev,nlat,nlon/), typeof(TWday35d))
   Gday0d_diff_high = new((/nlev,nlat,nlon/), typeof(Gday0d_sun))
   Gday0d_sun_high = new((/nlev,nlat,nlon/), typeof(Gday0d_sun))
   Gday0d_zero_high = new((/nlev,nlat,nlon/), typeof(Gday0d_sun))
   TWday35d_high = new((/nlev,nlat,nlon/), typeof(TWday35d))
   Gday0d_pop_sun_low = new((/nlev,nlat,nlon/), typeof(Gday0d_sun))
   Gday0d_pop_diff_low = new((/nlev,nlat,nlon/), typeof(Gday0d_sun))
   Gday0d_pop_zero_low = new((/nlev,nlat,nlon/), typeof(Gday0d_sun))
   TWday35d_pop_low = new((/nlev,nlat,nlon/), typeof(TWday35d))
   Gday0d_pop_diff_high = new((/nlev,nlat,nlon/), typeof(Gday0d_sun))
   Gday0d_pop_sun_high = new((/nlev,nlat,nlon/), typeof(Gday0d_sun))
   Gday0d_pop_zero_high = new((/nlev,nlat,nlon/), typeof(Gday0d_sun))
   TWday35d_pop_high = new((/nlev,nlat,nlon/), typeof(TWday35d))

   do i=0,nlon-1
    do j=0,nlat-1 
     do ii=0,nlev-1
     ii_25 = pctl_25(ii,j,i)
     ii_75 = pctl_75(ii,j,i)
     if(ismissing(ii_25) .or. ismissing(ii_75)) then; 
        continue
     end if
     ;print(ii_25)

     Gday_sun_low(ii,j,i) = Gday_sun_sort(ii_25,ii,j,i)
     Gday_sun_high(ii,j,i) = Gday_sun_sort(ii_75,ii,j,i)
     Gday_diff_low(ii,j,i) = Gday_diff_sort(ii_25,ii,j,i)
     Gday_diff_high(ii,j,i) = Gday_diff_sort(ii_75,ii,j,i)
     Gday_zero_low(ii,j,i) = Gday_zero_sort(ii_25,ii,j,i)
     Gday_zero_high(ii,j,i) = Gday_zero_sort(ii_75,ii,j,i)
     TWday_low(ii,j,i) = TWday_sort(ii_25,ii,j,i)
     TWday_high(ii,j,i) = TWday_sort(ii_75,ii,j,i)
     Gtwday_low(ii,j,i) = Gtwday_sort(ii_25,ii,j,i)
     Gtwday_high(ii,j,i) = Gtwday_sort(ii_75,ii,j,i)

     Gday0d_sun_low(ii,j,i) = Gday0d_sun_sort(ii_25,ii,j,i)
     Gday0d_sun_high(ii,j,i) = Gday0d_sun_sort(ii_75,ii,j,i)
     Gday0d_diff_low(ii,j,i) = Gday0d_diff_sort(ii_25,ii,j,i)
     Gday0d_diff_high(ii,j,i) = Gday0d_diff_sort(ii_75,ii,j,i)
     Gday0d_zero_low(ii,j,i) = Gday0d_zero_sort(ii_25,ii,j,i)
     Gday0d_zero_high(ii,j,i) = Gday0d_zero_sort(ii_75,ii,j,i)
     TWday35d_low(ii,j,i) = TWday35d_sort(ii_25,ii,j,i)
     TWday35d_high(ii,j,i) = TWday35d_sort(ii_75,ii,j,i)

     Gday0d_pop_sun_low(ii,j,i) = Gday0d_pop_sun_sort(ii_25,ii,j,i)
     Gday0d_pop_sun_high(ii,j,i) = Gday0d_pop_sun_sort(ii_75,ii,j,i)
     Gday0d_pop_diff_low(ii,j,i) = Gday0d_pop_diff_sort(ii_25,ii,j,i)
     Gday0d_pop_diff_high(ii,j,i) = Gday0d_pop_diff_sort(ii_75,ii,j,i)
     Gday0d_pop_zero_low(ii,j,i) = Gday0d_pop_zero_sort(ii_25,ii,j,i)
     Gday0d_pop_zero_high(ii,j,i) = Gday0d_pop_zero_sort(ii_75,ii,j,i)
     TWday35d_pop_low(ii,j,i) = TWday35d_pop_sort(ii_25,ii,j,i)
     TWday35d_pop_high(ii,j,i) = TWday35d_pop_sort(ii_75,ii,j,i)

     end do
    end do
   end do



   ;--note it's not medianingful to do ensemble statistics of area and population with Gday>0 at each grid cell 
   ;--rather it should use the ensemble statistics of Gday or TWday or WoE to determine whether a grid cell is impacted or not
    ;1) areal and population impact at each grid cell
    ;wgt3d_lev = conform_dims(dimsizes(Gday_sun_med),wgt,(/1,2/)) ; conform to 3d: [nlev, nlat, mlon]
    ;copy_VarCoords(Gday_sun_med, wgt3d_lev); land area is constant through time


    ;*focus all analysis on rural area and farmer population
    RA_lev = conform_dims(dimsizes(Gday_sun_med),rural_area,(/1,2/)) ; conform to 3d: [nlev, nlat, mlon]
    copy_VarCoords(Gday_sun_med, RA_lev); land area is constant through time

    ;***NOTE Outdoor sun scenario focuses on people who must work outdoors
    G0_area_sun_med := where(Gday_sun_med.gt.0, RA_lev, 0); grid land area with Gday_mx > 0 per year
    G0_pop_sun_med  := where(Gday_sun_med.gt.0, RP_lev_med, 0); focus on people who work on farms
    G0_area_sun_high := where(Gday_sun_high.gt.0, RA_lev, 0); grid land area with Gday_mx > 0 per year
    G0_pop_sun_high  := where(Gday_sun_high.gt.0, RP_lev_med, 0); 
    G0_area_sun_low := where(Gday_sun_low.gt.0, RA_lev, 0); grid land area with Gday_mx > 0 per year
    G0_pop_sun_low  := where(Gday_sun_low.gt.0, RP_lev_med, 0); 

    ;***NOTE to be comparable with outdoor sun, also mask shade/dark scenarios with rural population data
    G0_area_diff_med := where(Gday_diff_med.gt.0, RA_lev, 0); grid land area with Gday_mx > 0 per year
    G0_area_diff1_med := where(Gday_diff1_med.gt.0, RA_lev, 0); grid land area with Gday_mx > 0 per year
    G0_area_diff2_med := where(Gday_diff2_med.gt.0, RA_lev, 0); grid land area with Gday_mx > 0 per year
    G0_pop_diff_med  := where(Gday_diff_med.gt.0, RP_lev_med, 0); grid population with Gday_mx > 0 per year
    G0_pop_diff1_med  := where(Gday_diff1_med.gt.0, RP_lev_med, 0); grid population with Gday_mx > 0 per year
    G0_pop_diff2_med  := where(Gday_diff2_med.gt.0, RP_lev_med, 0); grid population with Gday_mx > 0 per year
    G0_area_zero_med := where(Gday_zero_med.gt.0, RA_lev, 0); grid land area with Gday_mx > 0 per year
    G0_pop_zero_med  := where(Gday_zero_med.gt.0, RP_lev_med, 0); grid population with Gday_mx > 0 per year
    TW35_area_med := where(TWday_med.gt.35, RA_lev, 0); grid land area with Gday_mx > 0 per year
    TW35_pop_med  := where(TWday_med.gt.35, RP_lev_med, 0); grid population with Gday_mx > 0 per year

    G0_area_diff_high := where(Gday_diff_high.gt.0, RA_lev, 0); grid land area with Gday_mx > 0 per year
    G0_pop_diff_high  := where(Gday_diff_high.gt.0, RP_lev_med, 0); grid population with Gday_mx > 0 per year
    G0_area_zero_high := where(Gday_zero_high.gt.0, RA_lev, 0); grid land area with Gday_mx > 0 per year
    G0_pop_zero_high  := where(Gday_zero_high.gt.0, RP_lev_med, 0); grid population with Gday_mx > 0 per year
    TW35_area_high := where(TWday_high.gt.35, RA_lev, 0); grid land area with Gday_mx > 0 per year
    TW35_pop_high  := where(TWday_high.gt.35, RP_lev_med, 0); grid population with Gday_mx > 0 per year

    G0_area_diff_low := where(Gday_diff_low.gt.0, RA_lev, 0); grid land area with Gday_mx > 0 per year
    G0_pop_diff_low  := where(Gday_diff_low.gt.0, RP_lev_med, 0); grid population with Gday_mx > 0 per year
    G0_area_zero_low := where(Gday_zero_low.gt.0, RA_lev, 0); grid land area with Gday_mx > 0 per year
    G0_pop_zero_low  := where(Gday_zero_low.gt.0, RP_lev_med, 0); grid population with Gday_mx > 0 per year
    TW35_area_low := where(TWday_low.gt.35, RA_lev, 0); grid land area with Gday_mx > 0 per year
    TW35_pop_low  := where(TWday_low.gt.35, RP_lev_med, 0); grid population with Gday_mx > 0 per year


  ;==summarize grid-cell specific areal and population impacts (median, low/high percentiles per warming level) to global 

   G0area_sun_med := tofloat(dim_sum_n(G0_area_sun_med, (/1,2/))*100*1e-6) ;sqrkm to Mha,total land area with Gday > 0 
   G0pop_sun_med := tofloat(dim_sum_n(G0_pop_sun_med, (/1,2/))*1e-6) ;million people
   G0area_diff_med := tofloat(dim_sum_n(G0_area_diff_med, (/1,2/))*100*1e-6) ;sqrkm to Mha,total land area with Gday > 0 
   G0area_diff1_med := tofloat(dim_sum_n(G0_area_diff1_med, (/1,2/))*100*1e-6) ;sqrkm to Mha,total land area with Gday > 0 
   G0area_diff2_med := tofloat(dim_sum_n(G0_area_diff2_med, (/1,2/))*100*1e-6) ;sqrkm to Mha,total land area with Gday > 0 
   G0pop_diff_med := tofloat(dim_sum_n(G0_pop_diff_med, (/1,2/))*1e-6) ;million people
   G0pop_diff1_med := tofloat(dim_sum_n(G0_pop_diff1_med, (/1,2/))*1e-6) ;million people
   G0pop_diff2_med := tofloat(dim_sum_n(G0_pop_diff2_med, (/1,2/))*1e-6) ;million people
   G0area_zero_med := tofloat(dim_sum_n(G0_area_zero_med, (/1,2/))*100*1e-6) ;sqrkm to Mha,total land area with Gday > 0 
   G0pop_zero_med := tofloat(dim_sum_n(G0_pop_zero_med, (/1,2/))*1e-6) ;million people
   TW35area_med := tofloat(dim_sum_n(TW35_area_med, (/1,2/))*100*1e-6) ;sqrkm to Mha,total land area with Gday > 0 
   TW35pop_med := tofloat(dim_sum_n(TW35_pop_med, (/1,2/))*1e-6) ;million people

   G0area_sun_low := tofloat(dim_sum_n(G0_area_sun_low, (/1,2/))*100*1e-6) ;sqrkm to Mha,total land area with Gday > 0
   G0pop_sun_low := tofloat(dim_sum_n(G0_pop_sun_low, (/1,2/))*1e-6) ;million people
   G0area_diff_low := tofloat(dim_sum_n(G0_area_diff_low, (/1,2/))*100*1e-6) ;sqrkm to Mha,total land area with Gday > 0
   G0pop_diff_low := tofloat(dim_sum_n(G0_pop_diff_low, (/1,2/))*1e-6) ;million people
   G0area_zero_low := tofloat(dim_sum_n(G0_area_zero_low, (/1,2/))*100*1e-6) ;sqrkm to Mha,total land area with Gday > 0
   G0pop_zero_low := tofloat(dim_sum_n(G0_pop_zero_low, (/1,2/))*1e-6) ;million people
   TW35area_low := tofloat(dim_sum_n(TW35_area_low, (/1,2/))*100*1e-6) ;sqrkm to Mha,total land area with Gday > 0
   TW35pop_low := tofloat(dim_sum_n(TW35_pop_low, (/1,2/))*1e-6) ;million people

   G0area_sun_high := tofloat(dim_sum_n(G0_area_sun_high, (/1,2/))*100*1e-6) ;sqrkm to Mha,total land area with Gday > 0
   G0pop_sun_high := tofloat(dim_sum_n(G0_pop_sun_high, (/1,2/))*1e-6) ;million people
   G0area_diff_high := tofloat(dim_sum_n(G0_area_diff_high, (/1,2/))*100*1e-6) ;sqrkm to Mha,total land area with Gday > 0
   G0pop_diff_high := tofloat(dim_sum_n(G0_pop_diff_high, (/1,2/))*1e-6) ;million people
   G0area_zero_high := tofloat(dim_sum_n(G0_area_zero_high, (/1,2/))*100*1e-6) ;sqrkm to Mha,total land area with Gday > 0
   G0pop_zero_high := tofloat(dim_sum_n(G0_pop_zero_high, (/1,2/))*1e-6) ;million people
   TW35area_high := tofloat(dim_sum_n(TW35_area_high, (/1,2/))*100*1e-6) ;sqrkm to Mha,total land area with Gday > 0
   TW35pop_high := tofloat(dim_sum_n(TW35_pop_high, (/1,2/))*1e-6) ;million people

   ;2) area-weighted average number of days per year with Gday>0 in impacted grid cells 
   ;only consider grid cells with Gday>0, set others to missing values
   wgt3d := conform_dims(dimsizes(Gday0d_sun_med),wgt,(/1,2/)) ; conform wgt to 3d: [nlev, nlat, mlon]
   copy_VarCoords(Gday0d_sun_med, wgt3d); gridcell area per warming level

   wgt_sun := where(Gday0d_sun_med .gt. 0, wgt3d, wgt3d@_FillValue)
   G0_avgdays_sun_med := tofloat(dim_sum_n(Gday0d_sun_med*wgt_sun, (/1,2/))/dim_sum_n(wgt_sun, (/1,2/)))
   G0_avgdays_sun_low := tofloat(dim_sum_n(Gday0d_sun_low*wgt_sun, (/1,2/))/dim_sum_n(wgt_sun, (/1,2/)))
   G0_avgdays_sun_high := tofloat(dim_sum_n(Gday0d_sun_high*wgt_sun, (/1,2/))/dim_sum_n(wgt_sun, (/1,2/)))
   wgt_diff := where(Gday0d_diff_med .gt. 0, wgt3d, wgt3d@_FillValue)
   G0_avgdays_diff_med := tofloat(dim_sum_n(Gday0d_diff_med*wgt_diff, (/1,2/))/dim_sum_n(wgt_diff, (/1,2/)))
   G0_avgdays_diff_low := tofloat(dim_sum_n(Gday0d_diff_low*wgt_diff, (/1,2/))/dim_sum_n(wgt_diff, (/1,2/)))
   G0_avgdays_diff_high := tofloat(dim_sum_n(Gday0d_diff_high*wgt_diff, (/1,2/))/dim_sum_n(wgt_diff, (/1,2/)))
   wgt_zero := where(Gday0d_zero_med .gt. 0, wgt3d, wgt3d@_FillValue)
   G0_avgdays_zero_med := tofloat(dim_sum_n(Gday0d_zero_med*wgt_zero, (/1,2/))/dim_sum_n(wgt_zero, (/1,2/)))
   G0_avgdays_zero_low := tofloat(dim_sum_n(Gday0d_zero_low*wgt_zero, (/1,2/))/dim_sum_n(wgt_zero, (/1,2/)))
   G0_avgdays_zero_high := tofloat(dim_sum_n(Gday0d_zero_high*wgt_zero, (/1,2/))/dim_sum_n(wgt_zero, (/1,2/)))
   wgt_TW := where(TWday35d_med .gt. 0, wgt3d, wgt3d@_FillValue)
   TW35_avgdays_med := tofloat(dim_sum_n(TWday35d_med*wgt_TW, (/1,2/))/dim_sum_n(wgt_TW, (/1,2/)))
   TW35_avgdays_low := tofloat(dim_sum_n(TWday35d_low*wgt_TW, (/1,2/))/dim_sum_n(wgt_TW, (/1,2/)))
   TW35_avgdays_high := tofloat(dim_sum_n(TWday35d_high*wgt_TW, (/1,2/))/dim_sum_n(wgt_TW, (/1,2/)))

   ;3) global heat stress burden (HSB): person days with Gday>0 
   G0_popdays_sun_med := tofloat(dim_sum_n(Gday0d_pop_sun_med, (/1,2/))*1e-6) ;million person days
   G0_popdays_diff_med := tofloat(dim_sum_n(Gday0d_pop_diff_med, (/1,2/))*1e-6) ;million person days
   G0_popdays_zero_med := tofloat(dim_sum_n(Gday0d_pop_zero_med, (/1,2/))*1e-6) ;million person days
   TW35_popdays_med := tofloat(dim_sum_n(TWday35d_pop_med, (/1,2/))*1e-6)
   G0_popdays_sun_low := tofloat(dim_sum_n(Gday0d_pop_sun_low, (/1,2/))*1e-6) ;million person days
   G0_popdays_diff_low := tofloat(dim_sum_n(Gday0d_pop_diff_low, (/1,2/))*1e-6) ;million person days
   G0_popdays_zero_low := tofloat(dim_sum_n(Gday0d_pop_zero_low, (/1,2/))*1e-6) ;million person days
   TW35_popdays_low := tofloat(dim_sum_n(TWday35d_pop_low, (/1,2/))*1e-6)
   G0_popdays_sun_high := tofloat(dim_sum_n(Gday0d_pop_sun_high, (/1,2/))*1e-6) ;million person days
   G0_popdays_diff_high := tofloat(dim_sum_n(Gday0d_pop_diff_high, (/1,2/))*1e-6) ;million person days
   G0_popdays_zero_high := tofloat(dim_sum_n(Gday0d_pop_zero_high, (/1,2/))*1e-6) ;million person days
   TW35_popdays_high := tofloat(dim_sum_n(TWday35d_pop_high, (/1,2/))*1e-6)

   G0_persons_sun_med  = tofloat(dim_sum_n(G0d_pop_sun_med, (/1,2/)))*1e-6 ;million person  
   G0_persons_diff_med  = tofloat(dim_sum_n(G0d_pop_diff_med, (/1,2/)))*1e-6 ;million person  
   G0_persons_zero_med  = tofloat(dim_sum_n(G0d_pop_zero_med, (/1,2/)))*1e-6 ;million person  
   G0_popavgdays_sun_med = G0_popdays_sun_med/G0_persons_sun_med
   G0_popavgdays_diff_med = G0_popdays_diff_med/G0_persons_diff_med
   G0_popavgdays_zero_med = G0_popdays_zero_med/G0_persons_zero_med

   ;after calculating ensemble median, there are also many numbers of 0<Gday0d<1, which after multiplying by population overestimates the impacted persons
   G0d_pop_sun_med := where(Gday0d_sun_med .ge. 1, G0d_pop_sun_med, G0d_pop_sun_med@_FillValue)
   G0d_pop_diff_med := where(Gday0d_diff_med .ge. 1, G0d_pop_diff_med, G0d_pop_diff_med@_FillValue)
   G0d_pop_zero_med := where(Gday0d_zero_med .ge. 1, G0d_pop_zero_med, G0d_pop_zero_med@_FillValue)
   Gday0d_pop_sun_med := where(Gday0d_sun_med .ge. 1, Gday0d_pop_sun_med, Gday0d_pop_sun_med@_FillValue)
   Gday0d_pop_diff_med := where(Gday0d_diff_med .ge. 1, Gday0d_pop_diff_med, Gday0d_pop_diff_med@_FillValue)
   Gday0d_pop_zero_med := where(Gday0d_zero_med .ge. 1, Gday0d_pop_zero_med, Gday0d_pop_zero_med@_FillValue)
   G0_persons_sun_med2  = tofloat(dim_sum_n(G0d_pop_sun_med, (/1,2/)))*1e-6 ;million person  
   G0_persons_diff_med2  = tofloat(dim_sum_n(G0d_pop_diff_med, (/1,2/)))*1e-6 ;million person  
   G0_persons_zero_med2  = tofloat(dim_sum_n(G0d_pop_zero_med, (/1,2/)))*1e-6 ;million person  
   G0_popdays_sun_med2  = tofloat(dim_sum_n(Gday0d_pop_sun_med, (/1,2/)))*1e-6 ;million person  
   G0_popdays_diff_med2  = tofloat(dim_sum_n(Gday0d_pop_diff_med, (/1,2/)))*1e-6 ;million person  
   G0_popdays_zero_med2  = tofloat(dim_sum_n(Gday0d_pop_zero_med, (/1,2/)))*1e-6 ;million person  
   G0_popavgdays_sun_med2 = G0_popdays_sun_med2/G0_persons_sun_med2
   G0_popavgdays_diff_med2 = G0_popdays_diff_med2/G0_persons_diff_med2
   G0_popavgdays_zero_med2 = G0_popdays_zero_med2/G0_persons_zero_med2

   id2 = get1Dindex(levels, (/1,2,4/)) 
   print("")
   print("total persons (farmer) with Gday>0 (sun) at "+levels(id2)+"C warming: "+G0pop_sun_med(id2)+" millions")
   print("total persons (farmer) with Gday0d>=1 (sun) at "+levels(id2)+"C warming: "+G0_persons_sun_med(id2)+" millions")
   print("total persons (farmer) with Gday0d>=1 (sun) at "+levels(id2)+"C warming: "+G0_persons_sun_med2(id2)+" millions")
   print("total person days (farmer) with Gday0d>=1 (sun) at "+levels(id2)+"C warming: "+G0_popdays_sun_med(id2)+" millions")
   print("total person days (farmer) with Gday0d>=1 (sun) at "+levels(id2)+"C warming: "+G0_popdays_sun_med2(id2)+" millions")
   print("population-weighted (farmer) average number of days with Gday0d>=1 (sun) at "+levels(id2)+"C warming: "+G0_popavgdays_sun_med2(id2))
   
   print("")
   print("total persons (farmer) with Gday>0 (shade) at "+levels(id2)+"C warming: "+G0pop_diff_med(id2)+" millions")
   print("total persons (farmer) with Gday0d>=1 (shade) at "+levels(id2)+"C warming: "+G0_persons_diff_med(id2)+" millions")
   print("total persons (farmer) with Gday0d>=1 (shade) at "+levels(id2)+"C warming: "+G0_persons_diff_med2(id2)+" millions")
   print("total person days (farmer) with Gday0d>=1 (shade) at "+levels(id2)+"C warming: "+G0_popdays_diff_med(id2)+" millions")
   print("total person days (farmer) with Gday0d>=1 (shade) at "+levels(id2)+"C warming: "+G0_popdays_diff_med2(id2)+" millions")
   print("population-weighted (farmer) average number of days with Gday0d>=1 (shade) at "+levels(id2)+"C warming: "+G0_popavgdays_diff_med2(id2))

   print("")
   print("total persons (farmer) with Gday>0 (dark) at "+levels(id2)+"C warming: "+G0pop_zero_med(id2)+" millions")
   print("total persons (farmer) with Gday0d>=1 (dark) at "+levels(id2)+"C warming: "+G0_persons_zero_med(id2)+" millions")
   print("total persons (farmer) with Gday0d>=1 (dark) at "+levels(id2)+"C warming: "+G0_persons_zero_med2(id2)+" millions")
   print("total person days (farmer) with Gday0d>=1 (dark) at "+levels(id2)+"C warming: "+G0_popdays_zero_med(id2)+" millions")
   print("total person days (farmer) with Gday0d>=1 (dark) at "+levels(id2)+"C warming: "+G0_popdays_zero_med2(id2)+" millions")
   print("population-weighted (farmer) average number of days with Gday0d>=1 (dark) at "+levels(id2)+"C warming: "+G0_popavgdays_zero_med2(id2))
 



  ;--------save ouput to netcdf and csv files


   lev1     =  levels
   lev1!0    = "level" ;named variables
   lev1@long_name = "global warming levels 1-4 degree Celsius"
  
   Gday_sun_med!0="lev"
   Gday_sun_med&lev = lev1
   Gday_sun_med!1="lat"
   Gday_sun_med&lat = lat
   Gday_sun_med!2="lon"
   Gday_sun_med&lon = lon
   copy_VarCoords(Gday_sun_med, Gday_diff_med)
   copy_VarCoords(Gday_sun_med, Gday_zero_med)
   copy_VarCoords(Gday_sun_med, Gday_sun_low)
   copy_VarCoords(Gday_sun_med, Gday_diff_low)
   copy_VarCoords(Gday_sun_med, Gday_zero_low)
   copy_VarCoords(Gday_sun_med, Gday_sun_high)
   copy_VarCoords(Gday_sun_med, Gday_diff_high)
   copy_VarCoords(Gday_sun_med, Gday_zero_high)
   copy_VarCoords(Gday_sun_med, TWday_med)
   copy_VarCoords(Gday_sun_med, TWday_low)
   copy_VarCoords(Gday_sun_med, TWday_high)
   copy_VarCoords(Gday_sun_med, RP_lev_med)
   copy_VarCoords(Gday_sun_med, RA_lev)

   copy_VarCoords(Gday_sun_med, G0_area_sun_med)
   copy_VarCoords(Gday_sun_med, G0_area_diff_med)
   copy_VarCoords(Gday_sun_med, G0_area_zero_med)
   copy_VarCoords(Gday_sun_med, G0_area_sun_low)
   copy_VarCoords(Gday_sun_med, G0_area_diff_low)
   copy_VarCoords(Gday_sun_med, G0_area_zero_low)
   copy_VarCoords(Gday_sun_med, G0_area_sun_high)
   copy_VarCoords(Gday_sun_med, G0_area_diff_high)
   copy_VarCoords(Gday_sun_med, G0_area_zero_high)
   copy_VarCoords(Gday_sun_med, TW35_area_med)
   copy_VarCoords(Gday_sun_med, TW35_area_low)
   copy_VarCoords(Gday_sun_med, TW35_area_high)
   copy_VarCoords(Gday_sun_med, G0_pop_sun_med)
   copy_VarCoords(Gday_sun_med, G0_pop_diff_med)
   copy_VarCoords(Gday_sun_med, G0_pop_zero_med)
   copy_VarCoords(Gday_sun_med, G0_pop_sun_low)
   copy_VarCoords(Gday_sun_med, G0_pop_diff_low)
   copy_VarCoords(Gday_sun_med, G0_pop_zero_low)
   copy_VarCoords(Gday_sun_med, G0_pop_sun_high)
   copy_VarCoords(Gday_sun_med, G0_pop_diff_high)
   copy_VarCoords(Gday_sun_med, G0_pop_zero_high)
   copy_VarCoords(Gday_sun_med, TW35_pop_med)
   copy_VarCoords(Gday_sun_med, TW35_pop_low)
   copy_VarCoords(Gday_sun_med, TW35_pop_high)

   copy_VarCoords(Gday_sun_med, Gday0d_sun_med)
   copy_VarCoords(Gday_sun_med, Gday0d_diff_med)
   copy_VarCoords(Gday_sun_med, Gday0d_zero_med)
   copy_VarCoords(Gday_sun_med, Gday0d_sun_low)
   copy_VarCoords(Gday_sun_med, Gday0d_diff_low)
   copy_VarCoords(Gday_sun_med, Gday0d_zero_low)
   copy_VarCoords(Gday_sun_med, Gday0d_sun_high)
   copy_VarCoords(Gday_sun_med, Gday0d_diff_high)
   copy_VarCoords(Gday_sun_med, Gday0d_zero_high)
   copy_VarCoords(Gday_sun_med, TWday35d_med)
   copy_VarCoords(Gday_sun_med, TWday35d_low)
   copy_VarCoords(Gday_sun_med, TWday35d_high)
   ;person-days with Gday>0 at each grid cell
   copy_VarCoords(Gday_sun_med, Gday0d_pop_sun_med)
   copy_VarCoords(Gday_sun_med, Gday0d_pop_diff_med)
   copy_VarCoords(Gday_sun_med, Gday0d_pop_zero_med)
   copy_VarCoords(Gday_sun_med, Gday0d_pop_sun_low)
   copy_VarCoords(Gday_sun_med, Gday0d_pop_diff_low)
   copy_VarCoords(Gday_sun_med, Gday0d_pop_zero_low)
   copy_VarCoords(Gday_sun_med, Gday0d_pop_sun_high)
   copy_VarCoords(Gday_sun_med, Gday0d_pop_diff_high)
   copy_VarCoords(Gday_sun_med, Gday0d_pop_zero_high)
   copy_VarCoords(Gday_sun_med, TWday35d_pop_med)
   copy_VarCoords(Gday_sun_med, TWday35d_pop_low)
   copy_VarCoords(Gday_sun_med, TWday35d_pop_high)


  ;===============save 3d data to netcdf
  dirNc  = dirPath2       ; where nc file will be saved
  filNc  = "Gday-gridcell-impacts-warminglevel-10yravg-ens-median-vcbc-outdoor-rural.nc"
  filNc2  = "rural_pop_warminglevels_12models-ens-median.nc" 

  pthNc  = dirNc+filNc
  system("/bin/rm -f "+pthNc)       ; remove any pre-existing file
  ncdf = addfile(pthNc ,"c")        ; open output netCDF file

  ncdf->Gday_sun_med  = Gday_sun_med
  ncdf->Gday_diff_med  = Gday_diff_med
  ncdf->Gday_zero_med  = Gday_zero_med
  ncdf->Gday_sun_low  = Gday_sun_low
  ncdf->Gday_diff_low  = Gday_diff_low
  ncdf->Gday_zero_low  = Gday_zero_low
  ncdf->Gday_sun_high  = Gday_sun_high
  ncdf->Gday_diff_high  = Gday_diff_high
  ncdf->Gday_zero_high  = Gday_zero_high
  ncdf->TWday_med  = TWday_med
  ncdf->TWday_low  = TWday_low
  ncdf->TWday_high  = TWday_high

  ncdf->G0_area_sun_med  = G0_area_sun_med
  ncdf->G0_area_diff_med  = G0_area_diff_med
  ncdf->G0_area_zero_med  = G0_area_zero_med
  ncdf->G0_area_sun_low  = G0_area_sun_low
  ncdf->G0_area_diff_low  = G0_area_diff_low
  ncdf->G0_area_zero_low  = G0_area_zero_low
  ncdf->G0_area_sun_high  = G0_area_sun_high
  ncdf->G0_area_diff_high  = G0_area_diff_high
  ncdf->G0_area_zero_high  = G0_area_zero_high
  ncdf->TW35_area_med  = TW35_area_med
  ncdf->TW35_area_low  = TW35_area_low
  ncdf->TW35_area_high  = TW35_area_high
  ncdf->G0_pop_sun_med  = G0_pop_sun_med
  ncdf->G0_pop_diff_med  = G0_pop_diff_med
  ncdf->G0_pop_zero_med  = G0_pop_zero_med
  ncdf->G0_pop_sun_low  = G0_pop_sun_low
  ncdf->G0_pop_diff_low  = G0_pop_diff_low
  ncdf->G0_pop_zero_low  = G0_pop_zero_low
  ncdf->G0_pop_sun_high  = G0_pop_sun_high
  ncdf->G0_pop_diff_high  = G0_pop_diff_high
  ncdf->G0_pop_zero_high  = G0_pop_zero_high
  ncdf->TW35_pop_med  = TW35_pop_med
  ncdf->TW35_pop_low  = TW35_pop_low
  ncdf->TW35_pop_high  = TW35_pop_high

  ncdf->Gday0d_sun_med  = Gday0d_sun_med
  ncdf->Gday0d_diff_med  = Gday0d_diff_med
  ncdf->Gday0d_zero_med  = Gday0d_zero_med
  ncdf->Gday0d_sun_low  = Gday0d_sun_low
  ncdf->Gday0d_diff_low  = Gday0d_diff_low
  ncdf->Gday0d_zero_low  = Gday0d_zero_low
  ncdf->Gday0d_sun_high  = Gday0d_sun_high
  ncdf->Gday0d_diff_high  = Gday0d_diff_high
  ncdf->Gday0d_zero_high  = Gday0d_zero_high
  ncdf->TWday35d_med  = TWday35d_med
  ncdf->TWday35d_low  = TWday35d_low
  ncdf->TWday35d_high  = TWday35d_high
  ncdf->Gday0d_pop_sun_med  = Gday0d_pop_sun_med
  ncdf->Gday0d_pop_diff_med  = Gday0d_pop_diff_med
  ncdf->Gday0d_pop_zero_med  = Gday0d_pop_zero_med
  ncdf->Gday0d_pop_sun_low  = Gday0d_pop_sun_low
  ncdf->Gday0d_pop_diff_low  = Gday0d_pop_diff_low
  ncdf->Gday0d_pop_zero_low  = Gday0d_pop_zero_low
  ncdf->Gday0d_pop_sun_high  = Gday0d_pop_sun_high
  ncdf->Gday0d_pop_diff_high  = Gday0d_pop_diff_high
  ncdf->Gday0d_pop_zero_high  = Gday0d_pop_zero_high
  ncdf->TWday35d_pop_med  = TWday35d_pop_med
  ncdf->TWday35d_pop_low  = TWday35d_pop_low
  ncdf->TWday35d_pop_high  = TWday35d_pop_high

  pthNc2  = dirNc+filNc2
  system("/bin/rm -f "+pthNc2)       ; remove any pre-existing file
  ncdf2 = addfile(pthNc2 ,"c")          ; open output netCDF file
  ncdf2->RP_lev_med = RP_lev_med; rural population 
  ncdf2->RA_lev_med = RA_lev; rural area

  ;===================Save 1d data to csv

  ;spatial statistics and extent of impact

  csv_file4 = dirPath2+"Gday0_area_pop_sun_warminglevels_10yravg_ens_median-vcbc-outdoor-rural.csv"
  system("rm -rf " + csv_file4)
  csv_file5 = dirPath2+"Gday0_area_pop_diff_warminglevels_10yravg_ens_median-vcbc-outdoor-rural.csv"
  system("rm -rf " + csv_file5)
  csv_file6 = dirPath2+"Gday0_area_pop_zero_warminglevels_10yravg_ens_median-vcbc-outdoor-rural.csv"
  system("rm -rf " + csv_file6)

  ;accumulative impact per year (days or person-days) after bias-correction of CMIP6 data
  csv_file7 = dirPath2+"Gday0_persdays_sun_warminglevels_10yravg_ens_median-vcbc-outdoor-rural.csv" 
  system("rm -rf " + csv_file7)
  csv_file8 = dirPath2+"Gday0_persdays_diff_warminglevels_10yravg_ens_median-vcbc-outdoor-rural.csv"
  system("rm -rf " + csv_file8)
  csv_file9 = dirPath2+"Gday0_persdays_zero_warminglevels_10yravg_ens_median-vcbc-outdoor-rural.csv"
  system("rm -rf " + csv_file9)

  csv_file11 = dirPath2+"TWday35_area_pop_warminglevels_10yravg_ens_median-vcbc-outdoor-rural.csv"
  system("rm -rf " + csv_file11)
  csv_file12 = dirPath2+"TWday35_persdays_warminglevels_10yravg_ens_median-vcbc-outdoor-rural.csv"
  system("rm -rf " + csv_file12)


;---Convert data array to 1D 
 
  ;2) extent of impact (area, population)
  field_names1 := (/"Level", "G0area_sun_med", "G0pop_sun_med","G0area_sun_high", "G0pop_sun_high","G0area_sun_low", "G0pop_sun_low"/) ;
  header1 := [/str_join(field_names1,",")/]
  write_table(csv_file4, "w", header1, "%s")
  alist1  := [/levels, G0area_sun_med,G0pop_sun_med,G0area_sun_high,G0pop_sun_high,G0area_sun_low,G0pop_sun_low/]
  format1 = "%g,%g,%g,%g,%g,%g,%g"
  write_table(csv_file4, "a", alist1, format1)

  ;3) average number of days with Gday>0
  field_names1 := (/"Level", "G0_avgdays_sun_med", "G0_popdays_sun_med","G0_avgdays_sun_high", "G0_popdays_sun_high", "G0_avgdays_sun_low", "G0_popdays_sun_low"/)
  header1 = [/str_join(field_names1,",")/]
  write_table(csv_file7, "w", header1, "%s")

  alist1  = [/levels, G0_avgdays_sun_med, G0_popdays_sun_med, G0_avgdays_sun_high, G0_popdays_sun_high, G0_avgdays_sun_low, G0_popdays_sun_low /]
  format1 = "%g,%g,%g,%g,%g,%g,%g"
  write_table(csv_file7, "a", alist1, format1)

  ;TW impacts
  field_names1 := (/"Level", "TW35area_med", "TW35pop_med","TW35area_high", "TW35pop_high","TW35area_low", "TW35pop_low"/)
  header1 := [/str_join(field_names1,",")/]
  write_table(csv_file11, "w", header1, "%s")
  alist1  := [/levels, TW35area_med,TW35pop_med, TW35area_high,TW35pop_high, TW35area_low,TW35pop_low/]
  format1 = "%g,%g,%g,%g,%g,%g,%g"
  write_table(csv_file11, "a", alist1, format1)

  field_names1 := (/"Level", "TW35_avgdays_med", "TW35_popdays_med", "TW35_avgdays_high", "TW35_popdays_high", "TW35_avgdays_low", "TW35_popdays_low"/) 
  header1 := [/str_join(field_names1,",")/]
  write_table(csv_file12, "w", header1, "%s")
  alist1  := [/levels, TW35_avgdays_med, TW35_popdays_med, TW35_avgdays_high, TW35_popdays_high, TW35_avgdays_low, TW35_popdays_low/]
  format1 = "%g,%g,%g,%g,%g,%g,%g"
  write_table(csv_file12, "a", alist1, format1)


  ;======diffuse (shade) scenario

  ;2) extent of impact (area, population)
  field_names1 := (/"Level", "G0area_diff_med", "G0pop_diff_med","G0area_diff_high", "G0pop_diff_high","G0area_diff_low", "G0pop_diff_low","G0area_diff1_med","G0pop_diff1_med","G0area_diff2_med","G0pop_diff2_med"/) ;
  header1 := [/str_join(field_names1,",")/]
  write_table(csv_file5, "w", header1, "%s")
  alist1  := [/levels, G0area_diff_med,G0pop_diff_med, G0area_diff_high,G0pop_diff_high, G0area_diff_low,G0pop_diff_low,G0area_diff1_med,G0pop_diff1_med,G0area_diff2_med,G0pop_diff2_med/]
  format1 = "%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g"
  write_table(csv_file5, "a", alist1, format1)

  ;3) annual number of days with Gday>0
  field_names1 := (/"Level", "G0_avgdays_diff_med", "G0_popdays_diff_med", "G0_avgdays_diff_high", "G0_popdays_diff_high", "G0_avgdays_diff_low", "G0_popdays_diff_low"/)
  header1 = [/str_join(field_names1,",")/]
  write_table(csv_file8, "w", header1, "%s")
  alist1  = [/levels, G0_avgdays_diff_med, G0_popdays_diff_med,G0_avgdays_diff_high, G0_popdays_diff_high,G0_avgdays_diff_low, G0_popdays_diff_low /]
  format1 = "%g,%g,%g,%g,%g,%g,%g"
  write_table(csv_file8, "a", alist1, format1)


  ;======zero radiation (dark) scenario

  ;2) extent of impact (area, population
  field_names1 := (/"Level", "G0area_zero_med", "G0pop_zero_med","G0area_zero_high", "G0pop_zero_high","G0area_zero_low", "G0pop_zero_low"/) ;
  header1 := [/str_join(field_names1,",")/]
  write_table(csv_file6, "w", header1, "%s")
  alist1  := [/levels, G0area_zero_med,G0pop_zero_med, G0area_zero_high,G0pop_zero_high, G0area_zero_low,G0pop_zero_low/]
  format1 = "%g,%g,%g,%g,%g,%g,%g"
  write_table(csv_file6, "a", alist1, format1)

  ;3) annual number of days with Gday>0
  field_names1 := (/"Level", "G0_avgdays_zero_med", "G0_popdays_zero_med", "G0_avgdays_zero_high", "G0_popdays_zero_high", "G0_avgdays_zero_low", "G0_popdays_zero_low"/)
  header1 = [/str_join(field_names1,",")/]
  write_table(csv_file9, "w", header1, "%s")
  alist1  = [/levels, G0_avgdays_zero_med, G0_popdays_zero_med,G0_avgdays_zero_high, G0_popdays_zero_high,G0_avgdays_zero_low, G0_popdays_zero_low /]
  format1 = "%g,%g,%g,%g,%g,%g,%g"
  write_table(csv_file9, "a", alist1, format1)



 end if; DATA






end 
